home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / cool / ge_cool.lha / GE_COOL2.1 / src / Matrix / Matrix.C < prev    next >
C/C++ Source or Header  |  1992-07-02  |  25KB  |  626 lines
  1. //
  2. // Copyright (C) 1991 Texas Instruments Incorporated.
  3. // Copyright (C) 1992 General Electric Company.
  4. //
  5. // Permission is granted to any individual or institution to use, copy, modify,
  6. // and distribute this software, provided that this complete copyright and
  7. // permission notice is maintained, intact, in all copies and supporting
  8. // documentation.
  9. //
  10. // Texas Instruments Incorporated, General Electric Company,
  11. // provides this software "as is" without express or implied warranty.
  12. //
  13.  
  14. #include <cool/Matrix.h>            // header of matrix and its envelope
  15.  
  16. // Matrix -- constructor specifiying size of Matrix
  17. // Input:    Row, Column parameters
  18. // Output:   None
  19.  
  20. template<class Type> 
  21. CoolMatrix<Type>::CoolMatrix<Type> (unsigned int rows, unsigned int cols)
  22. : CoolMatrix(rows, cols) {
  23.   this->data = new Type*[rows];            // Allocate the row memory
  24.   Type* columns = new Type[cols*rows];        // Allocate the array of elmts
  25.   for (int i = 0; i < rows; i ++)        // For each row in the Matrix
  26.     this->data[i] = &columns[i*cols];        // Fill in address of row
  27.   if (this->compare_s == NULL)            // If not set yet
  28.     this->compare_s = &Type##_CoolMatrix_is_data_equal; // Default is_equal
  29. }
  30.  
  31.  
  32. // Matrix -- constructor specifiying size of matrix and initial value
  33. // Input:    Row, Column parameters and initial value
  34. // Output:   None
  35.  
  36. template<class Type> 
  37. CoolMatrix<Type>::CoolMatrix<Type> (unsigned int rows, unsigned int cols,
  38.                     const Type& value)
  39. : CoolMatrix(rows, cols) {
  40.   this->data = new Type*[rows];            // Allocate the row memory
  41.   Type* columns = new Type[cols*rows];        // Allocate the array of elmts
  42.   for (int i = 0; i < rows; i ++) {        // For each row in the Matrix
  43.     this->data[i] = &columns[i*cols];        // Fill in address of row
  44.     for (int j = 0; j < cols; j++)        // For each element in column
  45.       this->data[i][j] = value;            // Assign initial value
  46.   }
  47.   if (this->compare_s == NULL)            // If not set yet
  48.     this->compare_s = &Type##_CoolMatrix_is_data_equal; // Default is_equal
  49. }
  50.  
  51. // Matrix -- constructor specifiying size of matrix and initial values
  52. // Input:    Rows, Column parameters and initial values
  53. // Output:   None
  54. // Note: Arguments in ... can only be pointers, primitive types like int, double,
  55. //       and NOT OBJECTS, passed by reference or value, like vectors, matrices;
  56. //       because constructors must be known and called at compile time!!!
  57. //       Cannot have char in ..., because char is 1 instead of 4 bytes, and 
  58. //       va_arg expects sizeof(Type) a multiple of 4 bytes.
  59.  
  60. template<class Type> 
  61. CoolMatrix<Type>::CoolMatrix<Type> (unsigned int rows, unsigned int cols, int n,
  62.                     Type v00 ...)
  63. : CoolMatrix(rows, cols) {
  64. #if ERROR_CHECKING
  65.   if (((sizeof(Type) % 4) != 0) ||        // Cause alignment problems
  66.       (sizeof(Type) > 8))            // User defined classes?
  67.     this->va_arg_error(#Type, sizeof(Type));    // So, cannot use this constructor
  68. #endif
  69.   this->data = new Type*[rows];            // Allocate the row memory
  70.   Type* columns = new Type[cols*rows];        // Allocate the array of elmts
  71.   for (int i = 0; i < rows; i ++)        // For each row in the Matrix
  72.     this->data[i] = &columns[i*cols];        // Fill in address of row
  73.   if (n > 0) {                    // If user specified values
  74.     va_list argp;                // Declare argument list
  75.     va_start (argp, v00);                   // Initialize macro
  76.     for (i = 0; i < rows && n; i++)        // For remaining values given
  77.       for (int j = 0; j < cols && n; j++, n--)    // Moving sequentially in Matrix
  78.     if ((i == 0) && (j == 0))
  79.       this->data[0][0] = v00;               // Hack for v00 ...
  80.     else
  81.       this->data[i][j] = va_arg(argp, Type); // Extract and assign
  82.     va_end(argp);
  83.   }
  84.   if (this->compare_s == NULL)            // If not set yet
  85.     this->compare_s = &Type##_CoolMatrix_is_data_equal; // Default is_equal
  86. }
  87.  
  88. // Matrix -- Construct a matrix from a block array of data, stored row-wise.
  89. // Input     number of rows and columns, and array of r*c data.
  90. // Ouput     None
  91.  
  92. template<class Type>
  93. CoolMatrix<Type>::CoolMatrix<Type> (unsigned int rows, unsigned int cols,
  94.                     const Type* data_block) 
  95. : CoolMatrix(rows,cols) {
  96.   this->data = new Type*[num_rows];        // Allocate the row memory
  97.   int n = num_rows * num_cols;            // Number of elements 
  98.   Type* columns = new Type[n];            // Allocate the array of elmts
  99.   for (int d = 0; d < n; d++)            // Copy all the data elements
  100.     columns[d] = data_block[d];
  101.   for (int i = 0; i < num_rows; i ++)         // For each row in the Matrix
  102.     this->data[i] = &columns[i*num_cols];    // Fill in address of row
  103. }
  104.  
  105.  
  106. // Matrix -- Copy constructor
  107. // Input:    Matrix reference
  108. // Output:   None
  109.  
  110. template<class Type> 
  111. CoolMatrix<Type>::CoolMatrix<Type> (const CoolMatrix<Type>& m)
  112. : CoolMatrix(m) {
  113.   this->data = new Type*[this->num_rows];    // Allocate the row memory
  114.   Type* columns = new Type[num_cols*num_rows];    // Allocate the array of elmts
  115.   for (int i = 0; i < num_rows; i ++) {        // For each row in the Matrix
  116.     this->data[i] = &columns[i*num_cols];    // Fill in address of row
  117.     for (int j = 0; j < this->num_cols; j++)    // For each element in column
  118.       this->data[i][j] = m.data[i][j];        // Copy value
  119.   }
  120. }
  121.  
  122.  
  123. // ~Matrix -- Destructor for Matrix class that frees up storage
  124. // Input:     *this
  125. // Output:    None
  126.  
  127. template<class Type> 
  128. CoolMatrix<Type>::~CoolMatrix<Type>() {
  129.   delete [] this->data[0];            // Free up the array of elmts
  130.   delete [] this->data;                // Free up the row memory
  131. }
  132.  
  133. // fill -- Set all elements of a matrix to a specified fill value
  134. // Input:  this*, reference to fill value
  135. // Output: None
  136.  
  137. template<class Type> 
  138. void CoolMatrix<Type>::fill (const Type& value) {
  139.   for (int i = 0; i < this->num_rows; i++)    // For each row in the Matrix
  140.     for (int j = 0; j < this->num_cols; j++)    // For each element in column
  141.       this->data[i][j] = value;            // Assign fill value
  142. }
  143.  
  144. // operator= -- Overload the assignment operator to assign a single 
  145. //              value to the elements of a Matrix. 
  146. // Input:       *this, reference to a value
  147. // Output:      Reference to updated Matrix object
  148.  
  149. template<class Type> 
  150. CoolMatrix<Type>& CoolMatrix<Type>::operator= (const Type& value) {
  151.   for (int i = 0; i < this->num_rows; i++)    // For each row in Matrix
  152.     for (int j = 0; j < this->num_cols; j++)    // For each column in Matrix
  153.       this->data[i][j] = value;            // Assign value
  154.   return *this;                    // Return Matrix reference
  155. }
  156.  
  157.  
  158. // operator= -- Overload the assignment operator to copy the elements
  159. //              in one Matrix to another. The existing storage for the 
  160. //              destination matrix is freed up and new storage of the same 
  161. //              size as the source is allocated.
  162. // Input:       *this, reference to Matrix
  163. // Output:      Reference to copied Matrix object
  164.  
  165. template<class Type> 
  166. CoolMatrix<Type>& CoolMatrix<Type>::operator= (const CoolMatrix<Type>& m) {
  167.   if (this != &m) {                // make sure *this != m
  168.     if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) {
  169.       delete [] this->data[0];            // Free up the array of elmts
  170.       delete [] this->data;            // Free up the row memory
  171.       this->num_rows = m.num_rows;        // Copy rows
  172.       this->num_cols = m.num_cols;        // Copy columns 
  173.       this->data = new Type*[this->num_rows];    // Allocate the rows
  174.       Type* columns = new Type[num_cols*num_rows]; // Allocate the columns
  175.       for (int i = 0; i < this->num_rows; i++)       // For each row
  176.     this->data[i] = &columns[i*num_cols];       // Fill in address of row
  177.     }
  178.     for (int i = 0; i < this->num_rows; i++)    // For each row in the Matrix
  179.       for (int j = 0; j < this->num_cols; j++)    // For each element in column
  180.     this->data[i][j] = m.data[i][j];    // Copy value
  181.   }
  182.   return *this;                    // Return Matrix reference
  183. }
  184.  
  185. // operator== -- Compare the elements of two Matrices of Type Type using
  186. //               the Compare pointer to funtion (default is ==). If one 
  187. //               Matrix has more rows and/or columns than the other, the
  188. //               result is FALSE
  189. // Input:        Reference to Matrix of Type Type
  190. // Output:       TRUE/FALSE
  191.  
  192. template<class Type> 
  193. Boolean CoolMatrix<Type>::operator== (const CoolMatrix<Type>& m) const {
  194.   if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
  195.     return FALSE;                            // Then not equal
  196.   for (int i = 0; i < this->num_rows; i++)                // For each row
  197.     for (int j = 0; j < this->num_cols; j++)                // For each columne
  198.       if ((*this->compare_s)(this->data[i][j],m.data[i][j]) == FALSE) // Same?
  199.     return FALSE;                              // Then no match
  200.   return TRUE;                              // Else same, so return TRUE
  201. }
  202.  
  203.  
  204. // is_data_equal -- Default data comparison function if user has not provided
  205. //                  another one. Note that this is not inline because we need
  206.     //                  to take the address of it for the compare static variable
  207. // Input:           Two Type references
  208. // Output:          TRUE/FALSE
  209.  
  210. template<class Type> CoolMatrix {
  211.   Boolean Type##_CoolMatrix_is_data_equal (const Type& t1, const Type& t2) {
  212.     return (t1 == t2);
  213.   }
  214. }
  215.  
  216.  
  217. // set_compare -- Specify the comparison function to be used
  218. //                in logical tests of vector elements
  219. // Input:         Pointer to a compare function
  220. // Output:        None
  221.  
  222. template<class Type> 
  223. void CoolMatrix<Type>::set_compare (Type##_CoolMatrix_Compare c) {
  224.   if (c == NULL)                // If no argument supplied
  225.     this->compare_s = &Type##_CoolMatrix_is_data_equal; // Default is_equal
  226.   else
  227.     this->compare_s = c;            // Else set to user function
  228. }
  229.  
  230. // operator<< -- Overload the output operator to print a matrix
  231. // Input:        ostream reference, Matrix reference
  232. // Output:       ostream reference
  233.  
  234. template<class Type> CoolMatrix {
  235.   ostream& operator<< (ostream& s, const CoolMatrix<Type>& m) {
  236.     for (int i = 0; i < m.rows(); i++) {    // For each row in matrix
  237.       for (int j = 0; j < m.columns(); j++)    // For each column in matrix
  238.     s << m.data[i][j] << " ";        // Output data element
  239.       s << "\n";                    // Output newline
  240.     }
  241.     return (s);                    // Return ostream reference
  242.   }}
  243.  
  244.  
  245.  
  246. // operator+= -- Destructive matrix addition of a scalar.
  247. // Input:        this*, scalar value
  248. // Output:       New matrix reference
  249.  
  250. template<class Type> 
  251. CoolMatrix<Type>& CoolMatrix<Type>::operator+= (const Type& value) {
  252.   for (int i = 0; i < this->num_rows; i++)    // For each row
  253.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  254.       this->data[i][j] += value;        // Add scalar
  255.   return *this;
  256. }
  257.  
  258. // operator*= -- Destructive matrix multiplication by a scalar.
  259. // Input:        this*, scalar value
  260. // Output:       New matrix reference
  261.  
  262.  
  263. template<class Type> 
  264. CoolMatrix<Type>& CoolMatrix<Type>::operator*= (const Type& value) {
  265.   for (int i = 0; i < this->num_rows; i++)    // For each row
  266.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  267.       this->data[i][j] *= value;        // Multiply by scalar
  268.   return *this;
  269. }
  270.  
  271. // operator/= -- Destructive matrix division by a scalar.
  272. // Input:        this*, scalar value
  273. // Output:       New matrix reference
  274.  
  275. template<class Type> 
  276. CoolMatrix<Type>& CoolMatrix<Type>::operator/= (const Type& value) {
  277.   for (int i = 0; i < this->num_rows; i++)    // For each row
  278.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  279.       this->data[i][j] /= value;        // division by scalar
  280.   return *this;
  281. }
  282.  
  283.  
  284. // operator+= -- Destructive matrix addition with assignment. Note that the
  285. //               dimensions of each matrix must be identical
  286. // Input:        this*, matrix reference
  287. // Output:       Updated this* matrix reference
  288.  
  289. template<class Type> 
  290. CoolMatrix<Type>& CoolMatrix<Type>::operator+= (const CoolMatrix<Type>& m) {
  291.   if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
  292.     this->dimension_error ("operator+=", #Type, 
  293.                this->num_rows, this->num_cols, m.num_rows, m.num_cols);
  294.   for (int i = 0; i < this->num_rows; i++)    // For each row
  295.     for (int j = 0; j < this->num_cols; j++)    // For each element in column
  296.       this->data[i][j] += m.data[i][j];        // Add elements
  297.   return *this;
  298. }
  299.  
  300.  
  301. // operator-= -- Destructive matrix subtraction with assignment. Note that the
  302. //               dimensions of each matrix must be identical
  303. // Input:        this*, matrix reference
  304. // Output:       Updated this* matrix reference
  305.  
  306. template<class Type> 
  307. CoolMatrix<Type>& CoolMatrix<Type>::operator-= (const CoolMatrix<Type>& m) {
  308.   if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
  309.     this->dimension_error ("operator-=", #Type, 
  310.                this->num_rows, this->num_cols, m.num_rows, m.num_cols);
  311.   int i, j;
  312.   for (i = 0; i < this->num_rows; i++)
  313.     for (j = 0; j < this->num_cols; j++)
  314.       this->data[i][j] -= m.data[i][j];
  315.   return *this;
  316. }
  317.  
  318. // operator* -- Non Destructive matrix multiply 
  319. //               num_cols of first matrix must match num_rows of second matrix.
  320. // Input:        two matrix references
  321. // Output:       New matrix containing the product.
  322.  
  323. template<class Type> CoolMatrix {
  324. CoolMatrix<Type>E operator* (const CoolMatrix<Type>& m1,
  325.                  const CoolMatrix<Type>& m2) {
  326.   if (m1.num_cols != m2.num_rows)        // dimensions do not match?
  327.     m1.dimension_error ("operator*=", #Type, 
  328.             m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  329.   CoolMatrix<Type> temp(m1.num_rows, m2.num_cols); // Temporary to store product
  330.   for (int i = 0; i < m1.num_rows; i++) {    // For each row
  331.     for (int j = 0; j < m2.num_cols; j++) {    // For each element in column
  332.       Type sum = 0;
  333.       for (int k = 0; k < m1.num_cols; k++)    // Loop over column values
  334.     sum += (m1.data[i][k] * m2.data[k][j]); // Multiply
  335.       temp(i,j) = sum;
  336.     }
  337.   }
  338.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  339.   return result;                     // copy of envelope
  340. }}
  341.  
  342. // operator- -- Non-destructive matrix negation. a = -b;
  343. // Input:       this*
  344. // Output:      New matrix 
  345.  
  346. template<class Type> 
  347. CoolMatrix<Type>E CoolMatrix<Type>::operator- () const {
  348.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  349.   int i, j;
  350.   for (i = 0; i < this->num_rows; i++)
  351.     for (j = 0; j < this->num_cols; j++)
  352.       temp.data[i][j] = - this->data[i][j];
  353.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  354.   return result;                     // copy of envelope
  355. }
  356.  
  357. // operator+ -- Non-destructive matrix addition of a scalar.
  358. // Input:       this*, scalar value
  359. // Output:      New matrix 
  360.  
  361. template<class Type> 
  362. CoolMatrix<Type>E CoolMatrix<Type>::operator+ (const Type& value) const {
  363.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  364.   for (int i = 0; i < this->num_rows; i++)    // For each row
  365.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  366.       temp.data[i][j] = (this->data[i][j] + value); // Add scalar
  367.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  368.   return result;                     // copy of envelope
  369. }
  370.  
  371.  
  372. // operator* -- Non-destructive matrix multiply by a scalar.
  373. // Input:       this*, scalar value
  374. // Output:      New matrix 
  375.  
  376. template<class Type> 
  377. CoolMatrix<Type>E CoolMatrix<Type>::operator* (const Type& value) const {
  378.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  379.   for (int i = 0; i < this->num_rows; i++)    // For each row
  380.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  381.       temp.data[i][j] = (this->data[i][j] * value); // Multiply
  382.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  383.   return result;                     // copy of envelope
  384. }
  385.  
  386.  
  387. // operator/ -- Non-destructive matrix division by a scalar.
  388. // Input:       this*, scalar value
  389. // Output:      New matrix 
  390.  
  391. template<class Type> 
  392. CoolMatrix<Type>E CoolMatrix<Type>::operator/ (const Type& value) const {
  393.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  394.   for (int i = 0; i < this->num_rows; i++)    // For each row
  395.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  396.       temp.data[i][j] = (this->data[i][j] / value); // Divide
  397.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  398.   return result;                     // copy of envelope
  399. }
  400.  
  401.  
  402. ////--------------------------- Additions ------------------------------------
  403.  
  404.  
  405. // transpose -- Return the transpose of this matrix.
  406. // Input:       this*
  407. // Ouput:       New matrix
  408.  
  409. template<class Type> 
  410. CoolMatrix<Type>E CoolMatrix<Type>::transpose() const {
  411.   CoolMatrix<Type> temp(this->num_cols, this->num_rows);
  412.   int i, j;
  413.   for (i = 0; i < this->num_cols; i++)
  414.     for (j = 0; j < this->num_rows; j++)
  415.       temp.data[i][j] = this->data[j][i];
  416.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  417.   return result;                     // copy of envelope
  418. }
  419.  
  420.  
  421. // abs -- Return the matrix of the absolute values.
  422. // Input:       this*
  423. // Ouput:       New matrix
  424.  
  425. template<class Type> 
  426. CoolMatrix<Type>E CoolMatrix<Type>::abs() const {
  427.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  428.   int i, j;
  429.   for (i = 0; i < this->num_rows; i++)
  430.     for (j = 0; j < this->num_cols; j++)
  431.       if (this->data[i][j] < 0)
  432.     temp.data[i][j] = - this->data[i][j];
  433.       else
  434.     temp.data[i][j] = this->data[i][j];
  435.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  436.   return result;                     // copy of envelope
  437. }
  438.  
  439. // sign -- Return the matrix whose elements are either -1,1 or 0
  440. // depending on whether the corresponding values are negative, positive, or 0.
  441. // Input:       this*
  442. // Ouput:       New matrix
  443.  
  444. template<class Type> 
  445. CoolMatrix<Type>E CoolMatrix<Type>::sign() const {
  446.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  447.   int i, j;
  448.   for (i = 0; i < this->num_rows; i++)
  449.     for (j = 0; j < this->num_cols; j++)
  450.       if (this->data[i][j] == 0)            // test fuzz equality to 0
  451.     temp.data[i][j] = 0;            // first.
  452.       else
  453.     if (this->data[i][j] < 0)
  454.       temp.data[i][j] = -1;
  455.     else
  456.       temp.data[i][j] = 1;
  457.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  458.   return result;                     // copy of envelope
  459. }
  460.  
  461. // element_product -- return the matrix whose elements are the products 
  462. // Input:     2 matrices m1, m2 by reference
  463. // Output:    New matrix, whose elements are m1[ij]*m2[ij].
  464.  
  465. template<class Type> CoolMatrix{
  466. CoolMatrix<Type>E element_product (const CoolMatrix<Type>& m1, const CoolMatrix<Type>& m2) {
  467.   if (m1.num_rows != m2.num_rows || m1.num_cols != m2.num_cols) // Size?
  468.     m1.dimension_error ("element_product", #Type, 
  469.             m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  470.   CoolMatrix<Type> temp(m1.num_rows, m1.num_cols);
  471.   int i, j;
  472.   for (i = 0; i < m1.num_rows; i++)
  473.     for (j = 0; j < m1.num_cols; j++)
  474.       temp.data[i][j] = m1.data[i][j] * m2.data[i][j];
  475.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  476.   return result;                     // copy of envelope
  477. }}
  478.  
  479. // element_quotient -- return the matrix whose elements are the quotients 
  480. // Input:     2 matrices m1, m2 by reference
  481. // Output:    New matrix, whose elements are m1[ij]/m2[ij].
  482.  
  483. template<class Type> CoolMatrix{
  484. CoolMatrix<Type>E element_quotient (const CoolMatrix<Type>& m1, const CoolMatrix<Type>& m2) {
  485.   if (m1.num_rows != m2.num_rows || m1.num_cols != m2.num_cols) // Size?
  486.     m1.dimension_error ("element_quotient", #Type, 
  487.             m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  488.   CoolMatrix<Type> temp(m1.num_rows, m1.num_cols);
  489.   int i, j;
  490.   for (i = 0; i < m1.num_rows; i++)
  491.     for (j = 0; j < m1.num_cols; j++)
  492.       temp.data[i][j] = m1.data[i][j] / m2.data[i][j];
  493.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  494.   return result;                     // copy of envelope
  495. }}
  496.  
  497. // update -- replace a submatrix of this, by the actual argument.
  498. // Input:       *this, starting corner specified by top and left.
  499. // Ouput:       mutated reference.
  500.  
  501. template<class Type> 
  502. CoolMatrix<Type>& CoolMatrix<Type>::update (const CoolMatrix<Type>& m, 
  503.                         unsigned int top, unsigned int left) {
  504.   unsigned int bottom = top + m.num_rows;
  505.   unsigned int right = left + m.num_cols;
  506.   if (this->num_rows < bottom || this->num_cols < right)
  507.     this->dimension_error ("update", #Type, 
  508.                bottom-top, right-left, m.num_rows, m.num_cols);
  509.   int i, j;
  510.   for (i = top; i < bottom; i++)
  511.     for (j = left; j < right; j++)
  512.       this->data[i][j] = m.data[i-top][j-left];
  513.   return *this;
  514. }
  515.  
  516.  
  517. // extract -- Return a submatrix specified by the top-left corner and size.
  518. // Input:       *this, starting corner specified by top and left, and size.
  519. // Ouput:       new matrix
  520.  
  521. template<class Type> 
  522. CoolMatrix<Type>E CoolMatrix<Type>::extract (unsigned int rows, unsigned int cols, unsigned int top, unsigned int left) const{
  523.   unsigned int bottom = top + rows;
  524.   unsigned int right = left + cols;
  525.   if ((this->num_rows < bottom) || (this->num_cols < right))
  526.     this->dimension_error ("extract", #Type, 
  527.                bottom-top, right-left, rows, cols);
  528.   CoolMatrix<Type> temp(rows, cols);
  529.   for (int i = 0; i < rows; i++)        // actual copy of all elements
  530.     for (int j = 0; j < cols ; j++)        // in submatrix
  531.       temp.data[i][j] = data[top+i][left+j];
  532.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  533.   return result;                     // copy of envelope
  534. }
  535.  
  536. // determinant -- Determinant of a square matrix using Cramer's rule.
  537. //              Signal Error exception if the matrix is not square.
  538.  
  539. template<class Type>
  540. Type CoolMatrix<Type>::determinant () const {
  541.   if (this->num_rows != this->num_cols || this->num_rows < 2)
  542.     this->dimension_error ("determinant", #Type,
  543.                this->num_rows, this->num_cols, 
  544.                this->num_rows, this->num_cols);
  545.   int n = this->num_rows, r, i, j;
  546.   Type det = 0, prod;
  547.   if (n == 2) {
  548.     det = (this->data[0][0] * this->data[1][1] - // border case of 2x2 matrix
  549.        this->data[0][1] * this->data[1][0]);
  550.   } else {
  551.     for (r = 0; r < n; r++) {            // compute sum of products
  552.       prod = 1;                    // along diagonals
  553.       for (i = r, j = 0; i < n; i++, j++)    // top-left to bottom-right
  554.     prod *= this->data[i][j];
  555.       for (i = 0; j < n; i++, j++)
  556.     prod *= this->data[i][j];
  557.       det += prod;                // coeft = +1
  558.     }
  559.     int e = n-1;                // index of last row/col
  560.     for (r = 0; r < n; r++) {            // compute sum of products
  561.       prod = 1;                    // along diagonals
  562.       for (i = r, j = e; i < n; i++, j--)    // top-right to bottom-left
  563.     prod *= this->data[i][j];
  564.       for (i=0; j >= 0; i++, j--)
  565.     prod *= this->data[i][j];
  566.     det -= prod;                // coeft = -1
  567.     }
  568.   }
  569.   return det;
  570. }
  571.  
  572. // dot_product -- Return the dot product of the row or column vectors
  573. // Input:       2 vectors by reference
  574. // Ouput:       dot product value
  575.  
  576. template<class Type> CoolMatrix {
  577.   Type dot_product (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  578.     if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols) // Size?
  579.       v1.dimension_error ("dot_product", #Type, 
  580.               v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  581.     Type dot = 0;
  582.     int i, j;
  583.     for (i = 0; i < v1.num_rows; i++)
  584.       for (j = 0; j < v1.num_cols; j++)        // generalized dot-product
  585.     dot += v1.data[i][j] * v2.data[i][j];    // of matrices
  586.     return dot;
  587. }}
  588.  
  589. // cross_2d -- Return the 2X1 cross-product of 2 2d-vectors
  590. // Input:       2 vectors by reference
  591. // Ouput:       cross product value
  592.  
  593. template<class Type> CoolMatrix {
  594. Type cross_2d (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  595.   if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols)
  596.     v1.dimension_error ("cross_2d", #Type, 
  597.             v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  598.   CoolMatrix<Type>& m1 = (CoolMatrix<Type>&) v1; // cast away const.
  599.   CoolMatrix<Type>& m2 = (CoolMatrix<Type>&) v2; 
  600.   return (m1.x() * m2.y()            // work for both col/row
  601.       -                    // representation.
  602.       m1.y() * m2.x());
  603. }}
  604.  
  605. // cross_3d -- Return the 3X1 cross-product of 2 3d-vectors
  606. // Input:       2 vectors by reference
  607. // Ouput:       3d cross product vector
  608.  
  609. template<class Type> CoolMatrix {
  610. CoolMatrix<Type>E cross_3d (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  611.   if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols)
  612.     v1.dimension_error ("cross_3d", #Type, 
  613.             v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  614.   CoolMatrix<Type> temp(v1.num_rows, v1.num_cols);
  615.   CoolMatrix<Type>& m1 = (CoolMatrix<Type>&) v1; // cast away const.
  616.   CoolMatrix<Type>& m2 = (CoolMatrix<Type>&) v2; 
  617.   temp.x() = m1.y() * m2.z() - m1.z() * m2.y(); // work for both col/row
  618.   temp.y() = m1.z() * m2.x() - m1.x() * m2.z(); // representation
  619.   temp.z() = m1.x() * m2.y() - m1.y() * m2.x();
  620.   CoolMatrix<Type>E& result = (CoolMatrix<Type>E&) temp; // same physical object
  621.   return result;                     // copy of envelope
  622. }}
  623.  
  624.  
  625.  
  626.